A  TUO-ENDED  SHOOTING  TECHNIQUE  FOR  CALCULATING  NORMAL 
NODES  IN  UNDERWATER  (U)  DEFENCE  RESEARCH  ESTABLISHMENT 
ATLANTIC  DARTMOUTH  (NOVA  SCOTI  D  D  ELLIS  SEP  85 
DREA-85/105  F/G  20/1  I 


AD-A161  965 


REPORT  85/105 
Scptcatwr  1985 


A  TWO-ENDED  SHOOTING  TECHNIQUE 
FOR  CALCULATING  NORMAL  MOOES 
IN  UNDERWATER  ACOUSTIC  PROPAGATION 

*  ;  -,V 

Dale  0.  Ellis 


Defence 

Research 

Establishment 

Atlantic 


Centre  de 
Recherches  pour  la 
Defense 
Atlantique 


CanadS 


UNLIMITED  DISTRIBUTION 


National  Defence 

Research  and 
Development  Branch 


Defense  Nationale 

Bureau  de  Recherche 
et  Developpment 


A  TWO-ENDED  SHOOTING  TECHNIQUE 
FOR  CALCULATING  NORMAL  MODES 
IN  UNDERWATER  ACOUSTIC  PROPAGATION 

Dale  D.  Ellis 


September  1985 


Approved  by  R.F.  Brown 

DISTRIBUTION  APPROVED  BY 


Director/Underwater  Acoustics  Division 

CHIEF  D.  R.E.A. 


REPORT  85/105 


Defence 

Research 

Establishment 

Atlantic 


Centre  de 
Recherches  pour  la 
Defense 
Atlantique 


Canad'a 


ABSTRACT 


An  algorithm  for  the  calculation  of  acoustic  normal  modes  in  the  ocean  is 
described.  The  algorithm  is  valid  for  an  arbitrary  sound  speed  and  density  profile 
in  the  water  column  and  bottom.  Losses  due  to  volume  absorption  in  the  water  and 
bottom,  surface  and  bottom  roughness,  and  shear  waves  in  the  bottommost  layer 
can  be  calculated  as  perturbations. 

The  essential  feature  of  the  algorithm  is  a  two-ended  shooting  method  in 
which  the  trial  solution  is  started  separately  at  the  surface  and  bottom  and 
numerically  integrated  to  a  matching  depth  in  the  middle,  usually  near  the 
minimum  sound  speed.  The  trial  solution  is  iterated  until  a  continuous  function 
with  a  continuous  derivative  is  obtained.  Shooting  from  both  ends  results  in  a  more 
stable  algorithm  and  gives  more  accurate  eigenfunctions  than  are  obtained  using 
conventional  single-ended  shooting  methods. 

This  paper  describes  the  theory  and  general  numerical  implementation  of 
the  algorithm.  For  completeness,  equations  are  given  for  calculating  group 
velocities,  propagation  loss  in  a  range-dependent  environment,  and  losses  due  to 
volume  absorption,  surface  and  bottom  roughness,  and  shear  waves  in  the 
bottommost  layer. 

The  two-ended  shooting  algorithm  has  been  implemented  in  the  computer 
program  PROLOS,  which  has  been  used  extensively  at  DREA  since  1979.  In  this 
paper  results  of  numerical  computations  will  be  presented,  including  the  AESD 
Workshop  Test  Cases  and  comparison  of  predictions  with  some  measured  data  from 
a  shallow  water  propagation  loss  experiment. 


RESUME 


On  decrit  un  algorithme  de  calcul  des  modes  acoustiques  normaux 
en  mer.  L*  algorithme  est  applicable  a  un  profil  arbitraire  de  density 
et  de  vitesse  du  son  dans  la  colonne  d'eau  et  au  fond.  Les 
affaiblissements  dus  a  l’absorption  par  la  masse  de  l'eau  et  du  fond,  k 
la  rugosity  de  la  surface  et  du  fond  et  aux  ondes  de  rotation  dans  la 
couche  inferieure  peuvent  etre  traites  comme  des  perturbations. 

L' element  essentiel  de  cet  algorithme  est  une  methode  de  "tir" 
bilateral  dans  laquelle  la  solution  d'essai  est  appliquee  separ6ment  k 
la  surface  et  au  fond  et  est  integree  numeriquement  jusqu'4  une 
profondeur  de  rencontre  au  milieu,  habituellement  pvka  du  point  de 
vitesse  du  son  minimum.  On  ameliore  par  iteration  la  solution  d'essai 
jusqu’a  ce  qu'on  obtienne  une  fonction  continue  avec  une  d6riv4e 
continue.  Le  "tir"  bilateral  donne  un  algorithme  plus  stable  et  des 
fonctions  propres  plus  precises  que  les  methodes  habituelles  de  "tir" 
unilateral. 

Le  present  document  contient  des  renseignements  sur  la  theorie  et 
le  calcul  numerique  general  de  1 ' algorithme .  Pour  plus  de  clarte,  on 
donne  les  equations  de  calcul  des  vitesses  de  groupe,  de 
l'affaiblissement  par  propagation  en  fonction  de  la  distance  et  des 
affaiblissements  dus  a  1' absorption  par  la  masse,  a  la  rugositd  de  la 
surface  et  du  fond  et  aux  ondes  de  rotation  dans  la  couche  inferieure. 

L' algorithme  de  "tir"  bilateral  a  fetfe  integre  au  programme 
informatique  PROLOS  qui  a  £t£  largement  utilise  au  CRDA  depuis  1979. 
Dans  le  present  document,  on  donne  les  resultats  des  calculs 
numeriques,  y  compris  ceux  des  cas  types  de  1' atelier  AESD,  et  on 
compare  des  valeurs  prdvues  avec  certaines  mesures  tir6es  d'une 
experience  d'affaiblissement  par  propagation  en  eau  peu  profonde. 
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1. 


INTRODUCTION 


To  predict  underwater  acoustic  propagation  from  a  knowledge  of  the 
environment,  a  number  of  models  -  primarily  ray  theoretic  and  normal  modes  - 
have  been  introduced.  Such  models  have  been  applied  to  a  variety  of  cases:  deep 
and  shallow  water,  high  and  low  frequency,  range-independent  and 
range-dependent  environments.  For  shallow  water  studies  normal  mode  models 
have  been  found  to  be  the  most  appropriate,  particularly  at  frequencies  below  « 
1  kHz.  Until  1978,  DREA  had  been  using  a  version  of  the  normal  mode  program  by 
Bartberger  and  Ackler  [1973]  of  the  United  States  Naval  Air  Development  Center. 
Since  this  program  was  written  originally  for  a  CDC  computer,  its  conversion  to 
run  on  a  DEC-20  system  revealed  two  numerical  problems  -  overflow  of  exponent 
and  loss  of  precision.  Both  problems  arose  from  the  shooting  technique  used  to 
solve  for  the  normal  mode  wave  numbers  and  mode  functions.  Although  the 
program  could  have  been  recoded  to  overcome  these  deficiencies,  an  improved 
algorithm  was  developed. 

This  paper  describes  the  alternate  algorithm  which  avoids  both  the  exponent 
overflow  and  the  loss  of  precision  of  the  Bartberger  shooting  algorithm  by  using  a 
two-ended  shooting  technique  instead  of  a  single-ended  method.  It  is  also  quite 
general  in  that  an  arbitrary  sound  speed  and  density  profile  can  be  handled.  Loss 
mechanisms,  such  as  volume  attenuation  in  the  water  and  bottom,  surface  and 
bottom  roughness,  and  shear  waves  can  be  incorporated  as  perturbations.  The 
formulation  is  given  in  terms  of  the  pressure  rather  than  the  velocity  potential, 
since  the  former  allows  density  changes  to  be  more  conveniently  handled,  and 
pressure  is  a  measurable  quantity. 

The  discussion  that  follows  will  first  give  a  short  description  of  normal  mode 
theory  and  present  the  equations  for  the  normal  modes  and  the  perturbative  losses. 
For  completeness  the  equations  for  calculating  group  velocities  and  propagation 
loss  in  a  range-dependent  environment  are  also  included.  This  is  followed  by  a 
description  of  the  two-ended  shooting  method  as  it  applies  to  normal  mode 
calculations  in  underwater  acoustics.  The  necessary  equations  are  presented  and 
cast  in  dimensionless  co-ordinates  for  use  in  computations. 

The  numerical  algorithm  can  be  implemented  in  a  number  of  ways.  Two 
computer  codes  have  been  developed  and  have  been  used  extensively  in  the  shallow 
water  studies  at  DREA.  One  [Ellis  and  Leverman,  1982]  uses  a  layered 
environmental  model  similar  to  Bartberger;  another  [MacEachem,  1983]  uses  a 
combination  of  the  Noumerov  method  and  a  Runge-Kutta  ordinary  differential 
equation  solver.  Results  of  numerical  computations  are  also  presented,  including 
the  AESD  Workshop  test  cases  [Spofford,  1973]  and  comparison  of  predictions  with 
some  measured  data  from  a  shallow  water  site. 
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2. 


THEORETICAL  BACKGROUND 


2.1  The  Normal  Mode  Equations 

The  environment  that  can  he  handled  by  normal  mode  theory  is  depicted  in 

Fig.  1.  It  consists  of  a  stratified  medium  in  which  the  sound  speed  and  density  can 

vary  arbitrarily  with  depth,  but  which  do  not  change  with  range.  In  its  exact  form 

the  boundaries  must  be  smooth  and  the  medium  must  not  be  lossy.  However,  it  is 

possible  using  perturbation  methods  to  extend  the  theory  to  include  the  effects  of 

volume  absorption  in  the  water  and  bottom,  scattering  from  rough  boundaries, 

weak  range-dependence,  and  shear  waves  in  the  bottommost  layer.  For  a  CW  point 

source  of  frequency  f  operating  at  range  r=0  and  depth  z  ,  the  pressure  field  P  at 

o 

an  arbitrary  depth  z  and  horizontal  range  r  from  the  source  is  given  by: 

P(r,z,t)  =  E  u  (z  )  u  (z)  H  (k  r)  e-lut  +  continuous  spectrum.  (1) 

n*l  n  0  n  on 

In  the  above  equation,  S  is  the  source  strength  (see  below),  p  is  the  density,  N  is  the 

number  of  trapped  modes,  the  u  are  the  normal  mode  functions,  H^  is  the  Hankel 

n  o 

function  of  zeroth  order  and  first  kind,  k  is  the  horizontal  wave  number  of  the  nth 

n 

mode,  o=2trf  is  the  angular  frequency  of  the  source,  and  t  is  the  time.  The  time 

dependence  is  not  important  for  what  follows  so  will  be  dropped  for  notational 

simplicity.  In  order  for  modes  to  be  trapped,  the  bottommost  layer  must  have  a 

sound  speed  which  is  greater  than  the  minimum  sound  speed.  A  rough  estimate  of 

the  number  of  trapped  modes  is  N=hf/c  ,  where  c  is  the  minimum  sound  speed 

m  zn 

and  h  is  the  water  depth.  A  better  approximation  [Ellis,  1982],  derived  using  WKB 
theory  [e.g.,  Schiff,  1968]  and  neglecting  density  changes,  is: 


z 


Fig.  1.  The  normal  mode  environmental  model.  Both  sound  speed  and 
density  can  vary  arbitrarily  and  discontinuously  with  depth. 
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N  =  [[V,  +  (2hfTi  /c  )(1  -  c  2/cJ)Vt ]] 
m  mu 


(2) 


where  [[x]]  is  used  to  denote  the  integer  of  part  of  x,  c  is  the  sound  speed  in  the 

B 

bottommost  layer  (extending  from  depth  z0  to  »),  and  tj  is  a  dimensionless  quantity 

B 

which  depends  only  on  the  shape  of  the  sound  speed  profile.  The  values  of  n  are 
usually  in  the  range  0.5  to  1.0.  As  examples,  t)=1  for  a  two-isoveloclty  layer 
model,  ti=0.79  for  a  parabolic  profile,  and  ri=0.67  for  a  bilinear  profile. 

If  N  >  1,  the  continuous  spectrum  usually  gives  negligible  contribution  to  the 
pressure  at  ranges  greater  than  several  water  depths,  and  is  usually  neglected  in 
calculations. 

There  are  a  number  of  conventions  for  what  is  meant  by  a  source  of  unit 
strength.  The  convention  implied  by  Eq.  (1)  is  that  the  pressure  is  unity  at  a  unit 
distance  from  the  source;  that  is,  near  the  source  the  pressure  is  given  by: 


P(x,t)  =  S|x-x  f 1  exp[i(k  |x-x  |  -  cat)]  (3) 

—  —  O  O  “  o 

where  x  =  (r,z),  k  =co/c  ,  and  c  is  the  sound  speed  at  the  source  position  x  .  S  has 
—  o  o  o  o 

numerical  value  of  unity,  but  has  the  appropriate  units  of  pressure  times  length. 

Once  the  pressure  has  been  obtained  using  Eq.  (1),  the  propagation  loss  can 
be  obtained  from: 

PL  =  -  10  l°gjP/Pref  I2  (4) 

where  P  is  the  pressure  at  •unit  distance  from  the  source.  Combining  Eqs.  (1),  (3) 

and  (4)  gives  the  expression  for  propagation  loss  in  terms  of  the  normal  modes  and 
wave  numbers: 

n  N  (1)  2 

PL(r,z)  =  -10  log  -  E  u  (z  )  u  (z)  H  (k  r)|  .  (5) 

10  0>  &.1  n  0  n  OD 

To  obtain  the  pressure  field,  Eq.  (1),  it  is  necessary  to  solve  the  normal 

2 

mode  equation  to  obtain  the  normal  mode  eigenvalues  k^  and  eigenfunctions  un<z). 
They  are  the  solutions  to  the  second  order  differential  equation: 

»  .  d  .  1  dll«.(z)  .  .  G)  .  2,  #  ,  .  t*\ 

p(z) — [ -  — S — ]  +  [ -  -  k  ]  u  (z)  »  0  (6) 

dz  p(z)  dz  c  2(z)  n  n 

and  satisfy  the  conditions: 
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(i)  u^(z)  is  continuous,  which  corresponds  physically  to  the  pressure  being 

continuous. 


(ii)  (du^/dz)/p(z)  is  continuous,  which  corresponds  to  the  vertical  particle 
velocity  being  continuous. 

(iii)  ^(01=0,  which  corresponds  to  zero  pressure  at  the  ocean  surface. 

(iv)  The  trapped  modes  have  finite  normalization.  If  the  bottommost  layer 
has  a  constant  density  and  sound  speed,  the  trapped  modes  behave  as: 

2  % 

u  (z)  ~  A  exp[(k  -  - )  z]  .  (7) 

n  n  c  a 

CB 


Under  the  above  conditions,  there  are  a  discrete  set  of  solutions  to  Eq.  (6), 

with  wave  numbers  k  constrained  by: 
n 


—  >  k  >  k  >  ...  >  k„  >  ii-  . 
c  12  N  c 

m  B 


(8) 


The  normal  mode  functions  are  orthogonal  with  respect  to  the  weighting  function 
p-1  and  are  normalized  to  unity,  that  is: 


J°°  — —  u  (z)  u  (z)  dz  =  6 
o  p(z)  m  n  mn 


(9) 


where  6  is  the  Krone  cker  delta. 

mn 


2.2  Modal  Loss  Calculations 

If  any  form  of  loss  is  introduced  into  the  system,  the  normal  mode  wave 
numbers  become  complex  with  a  small  imaginary  component;  that  is: 

K  *  k  +16  (10) 

n  n  n 

where  is  now  the  complex  wave  number.  The  expressions  for  the  pressure  and 

propagation  loss  are  unchanged  except  that  K  now  replaces  k  .  If  the  losses  are 

n  n 

small,  they  can  be  treated  as  perturbations  of  the  lossless  normal  mode  solutions. 
The  imaginary  part  6^  is  the  sum  of  the  various  attenuation  mechanisms: 

.  .bottom  t  water  tscatt  e  shear 

6=6  +6  +6  +6 


+  ... 


(11) 


where  ^ottom  is  the  loss  due  to  volume  absorption  in  the  bottom  layers,  6W&t'r  is 
n  n 

scatt 

the  loss  due  to  volume  absorption  in  the  water  column,  6^  is  the  loss  from  the 

shear 

coherent  field  due  to  scattering  from  rough  surfaces,  and  6q  is  the  loss  due  to 
shear  waves  in  the  bottommost  layer.  Expressions  for  the  losses  are  given  below. 


2.2.1  Bottom  absorption 

The  volume  absorption  coefficient  a  for  the  bottom  layers  experimentally 
seems  to  be  linearly  dependent  on  frequency  [Hamilton,  1974]  and  is  often  given  in 
the  units  dB/length  at  1  kHz.  Thus,  at  a  given  frequency,  the  absorption 
coefficient  is: 

c  .  Q2) 

20  log  e 
10 

which  has  units  of  nepers/length,  or  simply  inverse  length,  the  same  as  the  wave 

number  k  . 
n 

For  a  bottom  with  an  arbitrary  sound  speed,  density  and  absorption  profile, 
the  attenuation  coefficient  due  to  volume  absorption  is  given  by: 

.  ..  e(z)u2(z) 

6  bottom  =  (u/k  }  - 5 —  dz.  (13) 

n  n  h  p(z)  c(z) 

The  derivation  of  this  formula  is  given  in  Appendix  B.  The  attenuation  due  to  any 
layer  can  be  obtained  by  restricting  the  range  of  integration  to  that  layer.  If  the 
bottommost  layer  supports  shear  waves,  then  it  is  not  included  in  the  integration; 
the  equations  of  Sec.  2.2.4  are  used  instead. 


2.2.2  Volume  absorption  in  the  water  column 

The  attenuation  due  to  absorption  of  acoustic  energy  by  the  water  column 
can  be  treated  in  an  identical  manner  to  the  absorption  by  the  bottom  layers. 

However,  the  plane  wave  absorption  coefficient  varies  as  f2  rather  than  linearly  as 
for  the  sediments.  The  absorption  coefficient  at  any  frequency  can  be  obtained 
from  a  formula  such  as  Thorpe's  [1967]: 


r  0-1  f  .  40  f  .  10 

a  *  I - +  - J  - 

W  1  +  (f/1000)  2  4100  +  (f/1000)  2  0.9144 


has  units  of  dB/meter.  The  corresponding  modal  attenuation  coefficient  is  then 
(Appendix  B): 


h  \ 

I  H 


.u  2(z) 

6  water  _  (a  £  f  k  ,  jh  _n -  dz 

n  w  w  n  o  c(2) 


where  p^  is  the  density  of  water  and  where: 


c  =  a  /20  log  e. 

W  W  10 


(15) 


2.2.3  Surface  and  bottom  scattering 

Scattering  losses  due  to  rough  boundaries  can  be  included  in  an  approximate 
way.  By  ignoring  the  contribitions  from  the  non- specularly  reflected  energy,  and 
using  the  Kirchhoff  approximation,  Kuperman  and  Ingenito  [1977]  find  the  modal 
attenuation  coefficient  for  surface  and  bottom  roughness  to  be: 


^S-scatt 

n 


go  Vo) 

2p(o)kQ 


rdun{°) 


L  dz  J 


^B-scatt 

°n 


<  V*' 

2p(h~  )k 


du  (h  ) 
n 


dz 


♦  lY. 


(16) 


(17) 


where 


V2>  - 


(g>  2/c  2(z)  -  k^) 


V. 


a  and  a,  are  the  r.m.s.  surface  and  bottom  roughness  heights,  and  the 

O  1 

notation  h  means  that  the  various  functions  are  to  be  evaluated  Just  above  the 

water-bottom  interface.  Also  note  that  if  o/c(o)  5  k  or  «/c(h  )  5  k  ,  then  the 

n  n 

corresponding  scattering  losses  are  zero. 


These  formulae  are  valid  in  the  small  wave  height  approximation: 


2  a2  y3  (o)  «  1  , 
o  'n 

and  2  c  2  y2  (h  )  «  1. 
i  n 


(18) 
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2.2.4  Shear  waves  in  the  bottommost  layer 

When  the  speed  of  shear  waves  in  the  bottommost  layer  c  is  considerably 

s 

less  than  the  minimum  speed  of  the  compressional  waves  at  all  depths,  the  losses 
can  be  treated  by  perturbation  methods.  A  formula  is  given  in  a  report  on  the 
SNAP  normal  mode  program  [Jensen  and  Ferla,  1979]: 


^  shear 
n 


2,  +  \ 
"n  (ZB  1 


0 

i 

"8k" 


U  ♦ 


/  \ \ 2  Q 

—  JL  ]  qo  ) 

P  0  1 


(19) 


where  z*  refer  to  points  Just  above  (-)  and  Just  below  (+)  the  top  of  the  bottommost 
layer, 


0 


2 

O 


c3(V 


» 


0 


2 


G)2  )  % 

c  a(z  J) 


and 


Q(0  )  = 

i 


1-IRI2 


2  Imag(R) 


if  Real  0^a  i  0 
if  Real  0X  2  <  0  . 


(20) 


R(B^)  is  the  plane  wave  reflection  coefficient  [Brekhovskikh,  1980], 


R(0  )  * 

i 


p  0  C  +  p  (0  0  /0  )S 

2  1  2S  2  12  S  2S 


p  0  C  +p(0  0/0)S  +P0 

2  1  2S  2  12  S  25  i: 


(21) 


where 


» 


7 


C  -  [1  -  2(k  c  /«)  ] 
2s  ns 


S 

25 


-  c 


2S' 


Note  that  this  formula  contains  losses  due  to  both  the  compressional  and 
shear  wave  losses. 


2.3  Group  Velocities 

The  speed  at  which  acoustic  energy  is  transported  in  mode  n  at  frequency  a 
is  given  by  the  group  velocity  v  =  |du/dk  |.  This  expression  is  often  evaluated  by 

o  ^ 

solving  the  normal  mode  Eq.  (6)  at  two  frequencies  o±Ao  and  forming  the  finite 

difference  estimate  v  u  =  |Aco/Ak  |.  However,  in  addition  to  the  extra 

g  n 

computational  effort,  the  technique  is  subject  to  numerical  error. 

An  exact  and  alternate  expression  for  the  group  velocity  [Tolstoy,  1956; 
Koch  et  al,  1983]  is  given  by: 


v 


(n) 

g 


u2(z) 

n 

p(z)  c2  (z) 


(22) 


This  formula  is  similar  to  the  expression  for  the  volume  attenuation  coefficient 
[Appendix  B],  and  can  be  derived  in  a  similar  manner  [Chapman  and  Ellis,  1983]. 
Because  of  the  similarity  of  the  above  integral  with  Eqs.  (13)  and  (15),  the  group 
velocity  can  be  obtained  with  essentially  no  additional  computation. 


2.4  Extensions  for  Comparison  with  Data 

The  expressions  for  the  pressure  and  propagation  loss,  Eqs.  (1)  and  (5),  can 
be  easily  extended  to  facilitate  comparison  with  measurements  in  a  weakly 
range -dependent  environment  and  for  a  broadband  source. 

The  expression  for  the  propagation  loss  given  in  Eq.  (5)  is  a  coherent 
summation  of  the  modes  that  includes  all  the  interference  terms.  It  is  sometimes 
more  appropriate  to  sum  the  modes  incoherently  to  obtain: 

N  m 

IPL(r.z)  =  -  10  log  {[tt/p(z  )J 2  E  |u  (z  )  u  (z)  Hu,(k  r)|  2 }  . 

xo  o  .non  on 
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This  is  a  smooth  function  of  range  since: 


|H(1)(k  r)|  2  ~  (2/irk  r)  e"26  nr 
on  n 


where  the  asymptotic  form  of  the  Hanke  function  has  been  used,  and  the  mode 

attenuation  6  ,  Eq.  (10),  has  been  included.  The  neglected  terms  in  the  propagation 
n 

loss  are  of  the  form: 


N 

(2/irr)  £ 

n,m=l 
n*m 


u  (z  )  u  (2)  u(z)  u  (z)  Ik  kj 
non  mom  nm 


V, 


exp[i(k  -k  )r  -  (6  +6 Jr] 
n  m  n  m 


which  are  expected  to  be  small  for  a  number  of  reasons.  First,  the  oscillating 

terms  exp[i(k  -k  )r]  will  tend  to  make  the  range-averaged  value  small.  Second,  if 
n  m 

averaged  over  a  frequency  band,  there  will  be  additional  smoothing. 

The  incoherent  propagation  loss  should  be  used  cautiously  for  deep  water 
environments  since  the  regions  of  high  and  low  intensity  including  convergence 
zones  are  averaged  out.  However,  in  shallow  water  it  is  a  useful  quantity,  and  for 
broadband  sources  is  more  meaningful  than  the  coherent  propagation  loss. 

For  an  environment  in  which  the  sound  speed  profile  and  water  depth  are 
slowly  varying,  mode  coupling  can  be  neglected  and  the  adiabatic  approximation 
holds  [Pierce,  1965].  The  form  of  the  expressions  for  the  pressure  and  propagation 
loss  is  unchanged  except  that: 

(i)  u  (z  )  is  obtained  from  the  environment  at  range  r*0; 

n  o 

(ii)  u  (z)  is  obtained  from  the  environment  at  range  rs 

n 

(iii)  k  is  the  average  value  of  the  wave  number  in  the  range  0  to  r;  and, 

n 

(iv)  the  number  of  modes  is  determined  from  the  minimum  value  of  N  at  all 
ranges  between  0  and  r. 

In  principle,  the  normal  modes  should  be  calculated  at  every  range  for  which  the 
propagation  loss  is  desired.  In  practice,  the  mode  functions  and  wave  numbers  are 
evaluated  at  a  relatively  small  number  of  ranges  and  interpolated  for  intermediate 
ranges. 


3.  NUMERICAL  IMPLEMENTATION 

A  common  procedure  for  solving  for  the  normal  mode  eigenfunctions  and 

eigenvalues  is  to  numerically  integrate  the  normal  mode  Eq.  (6)  using  a  shooting 

method.  In  this  method  a  trial  value  is  chosen  for  k  ,  then  starting  at  some  large 

n 

value  of  z  deep  in  the  ocean  bottom  with  suitable  values  of  u  and  du/dz,  the 

solution  is  numerically  integrated  toward  the  ocean  surface41  at  z=0.  In  general, 

the  solution  obtained  from  the  trial  value  of  k  will  not  satisfy  the  condition 

n 

u  (0)=0,  but  based  on  the  calculated  values  of  u  and  u'  at  z=0  a  new  value  of  k  can 
n  n 

be  estimated.  The  procedure  is  repeated  until  the  value  of  u  at  z=0  is  sufficiently 

n 

small,  or  until  k  has  the  required  degree  of  accuracy. 


3.1  Some  Difficulties  With  Shooting  Methods 

Two  numerical  problems  can  occur  with  this  method  of  shooting,  namely, 
exponent  overflow  and  loss  of  precision.  In  solving  the  normal  mode  equation  in 
underwater  acoustics,  exponent  overflow  will  occur  if  the  starting  point  is  chosen 
too  deep,  that  is  if  z  is  too  large.  Although  less  likely  near  the  surface,  it  can  also 

occur  there  if  the  sound  speed  increases  near  the  surface,  that  is  if  c>o/k  .  The 

n 

increase  in  sound  speed  near  the  surface  can  cause  a  second  problem  as  well  -  loss 
of  precision  -  because  the  true  solution  gets  increasingly  buried  in  the  roundoff  and 
truncation  errors  as  the  integration  proceeds  toward  the  ocean  surface.  These 
difficulties  will  be  discussed  later  in  more  detail,  and  a  method  will  be  proposed  for 
getting  around  these  problems. 


Fig.  2.  Normal  mode  solution  u  (z)  (dashed  line)  superimposed  on  a 
sound  speed  profile  c(z)  (solid  line). 


♦One  cannot  start  at  the  surface  and  integrate  downward  into  the  bottom  because 
numerical  roundoff  would  cause  the  functions  to  grow  exponential ’y  even  for  the 
"exact"  k  where  the  functions  should  decay  exponentially. 
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Before  elaborating  on  these  points,  however,  it  is  useful  to  look  at  some 
properties  of  the  mode  function.  Fig.  2  shows  a  mode  function  superimposed  on  a 

sound  speed  profile  at  the  phase  velocity  v  =  u/k  of  the  mode.  A  number  of 

n  n 

properties  of  the  mode  function  to  note  are  as  follows:  from  the  boundary  condition 
at  the  surface, 

(1)  u  (0)=0; 

n 

and  from  y  as  defined  following  Eq.  (17), 

H 

(2)  for  2  <  Zj,  then  <  0>  and  u^z)  has  an  increasing  exponential  type  of 
behaviour.  The  solution  is  approximately  given  by: 


where. 


u  (z)  =  exp[I  (z)]  -  exp[-I  (z)] 
n  i  i 


Hz)  =  J0hrn(z)|  dz 


(3)  for  z^<  z  <  z^,  then  u^(z)  is  oscillatory.  The  solution  is  approximately 


given  by: 


where, 


u  (z)  *  A(z)  sind  (z)  +  <J>) 
n  a 


Hz)  -  I  Yn(z)dz 
x 

A(z)  is  slowly  varying,  and  <J>  is  the  phase,  approximately  ir/4,  at  the 
turning  point  z  . 

(4)  for  z  >  z  ,  u  (z)  has  a  decreasing  exponential  type  of  behaviour.  The 
2  n 

boundary  condition  at  infinity  requires  that  only  the  decreasing 
exponential  be  present;  thus, 

u  (z)  .  C  exp  [-1  <z)]  (25) 

n  a 

where  C  is  a  real  constant,  and 


Uz)  *  Jz  |y^(z)|  dz. 


(5)  In  the  interval  z  <  z  <  z  ,  u  (z)  has  n-1  zeros, 
i  an 


-  ^ 
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The  numerical  problems  occur  in  the  exponential  regions.  Suppose  we 
started  at  z=0  and  integrated  the  differential  equation  downwards  toward  z  . 

Although  the  increasing  and  decreasing  terms  have  equal  amplitude  at  z=0,  the 
increasing  exponential  soon  dominates  and  the  second  term  gets  lost  in  the 
roundoff  error.  Any  error  that  occurs  in  the  coefficient  of  exp[-I  (z)]  is  soon 

damped  out,  and  the  method  is  stable.  The  only  problem  that  might  occur  is  that 
I(z)  becomes  large  enough  so  that  exp[I(z)]  may  cause  overflow  in  computer 
calculations.  On  the  other  hand,  suppose  that  we  have  the  correct  solution  at  z*z^ 

and  want  to  integrate  (shoot)  upwards  toward  z=0;  the  procedure  leads  to  loss  of 
precision.  Suppose  there  is  a  slight  roundoff  so  that  the  coefficient  of  exp[-Mz)]  is 

no  longer  unity  but  1+c,  where  c  is  a  very  small  quantity.  Even  with  no  further 
roundoff,  the  solution  is  no  longer  zero  but: 

u  (o)  =  -c  exp[I  (z  )]  (27) 

n  11 

so  we  have  lost  I  (z  )ln(10)  significant  digits  in  the  accuracy  of  the  mode 

function.  This  in  turn  will  cause  problems  in  the  algorithm  that  predicts  the  next 

trial  value  of  k  .  [The  final  accuracy  of  k  is  not  necessarily  degraded,  since  the 
n  n 

numerical  procedure  can  be  terminated  as  long  as  k  is  sufficiently  accurate,  even 

n 

though  the  condition  u  (0)  <  tolerance  is  not  satisfied.  However,  convergence  may 

n 

be  slower,  and  the  mode  function  will  be  less  accurate  for  z  <  z  ;  that  is,  near  the 

i 

surface  where  solutions  are  very  often  needed.]  This  loss  of  accuracy  can  be 

eliminated  by  integrating  downward  from  z=0.  Similarly,  the  solution  for  large 

values  of  z  can  be  accurately  obtained  by  integrating  upward  toward  the  surface. 

The  numerical  problems  do  not  occur  for  z  >  z  since  we  are  integrating  upwards 

2 

and  any  part  of  the  unwanted  solution  gets  quickly  damped  out. 

The  method  of  matching  the  upward  and  downward  solutions  and  getting  an 
improved  estimate  for  the  eigenvalue  is  discussed  in  the  next  section. 


3.2  The  Two-Ended  Shooting  Method 

The  previously  discussed  loss  of  precision  encountered  in  the  single-ended 
shooting  method  can  be  remedied  by  using  a  two-ended  shooting  method.  The 
problem  of  exponent  overflow  in  numerical  computations  is  easily  handled,  so 
discussion  is  postponed  until  later.  We  first  discuss  the  two-ended  shooting  method 
and  the  iteration  procedure  for  refining  the  accuracy  of  the  trial  eigenfunction,  or 
wave  number. 

The  method  described  here  is  a  generalization  of  that  of  Blatt  [1967]  for  the 
solution  of  the  Schrodinger  equation  of  quantum  mechanics.  The  generalization 
allows  for  density  changes  in  the  acoustic  medium,  a  feature  that  does  not  occur  in 
the  quantum  mechanical  problem.  The  generalized  algorithm  which  is  used  to 
refine  the  approximate  eigenvalue  is  obtained  from  an  application  of  the 
variational  principle,  and  is  derived  in  Appendix  C. 
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3.2.1  Outline  of  the  shooting  method 


A  brief  sketch  of  the  two- ended  shooting  method  will  now  be  presented 
using  Fig.  2  for  reference.  First,  a  matching  depth  is  chosen  between  the 

turning  points  z  and  z  ;  a  convenient  choice  for  this  depth  is  the  minimum  sound 
12 

speed,  since  the  solution  is  always  oscillatory  at  the  minimum  sound  speed.  Next, 
an  estimate  is  made  for  the  eigenvalue  k^  and  the  turning  points  z^  and  z^  are 
determined.  Then,  suitable  limits  of  integration  z  and  z.  are  determined.  In 

ft  D 

practice  these  limits  are  often  the  ocean  surface  and  the  top  of  the  bottommost 
layer,  respectively.  Suitable  boundary  conditions  are  chosen  at  z  and  z^  and  the 

differential  equation  is  integrated  (down)  from  z  to  the  matching  depth  z  ,  and 

a  m 

(up)  from  z.  to  z  .  The  integration  procedure  gives  the  function  u  and  its 

D  !Q  & 

derivative  at  z  from  the  two  integrations.  Since  the  normalizations  are 
m 

arbitrary,  one  of  the  mode  functions  can  be  renormalized  to  give  a  function  which 
is  continuous,  but  which  in  general  has  a  discontinuous  derivative.  Provided  the 
number  of  zero  crossings  of  the  mode  function  is  correct,  an  improved 

approximation  for  the  eigenvalue  ka  is  given  by  (Appendix  C): 


k  2 
n 


k 1 « 


du  (z  ) 
n  m 

dz 


P<  *1) 

m 


du  (z  ) 
n  m 

dz 


u  (z  ) 
n  m 


z.  u  2 (z) 
b  n 


piz) 


dz 


(28) 


The  above  equation  is  the  key  formula  for  the  improvement  of  the  eigenvalue.  The 
notation  z~  or  z+  means  that  the  functions  p  or  u  are  to  be  evaluated  Just  above 

IQ  XZ1  XI 

(-)  or  Just  below  (+)  the  matching  depth  z  .  The  procedure  is  repeated  until  the 

m 

correction  term  is  small  enough,  and/or  until  the  derivatives  match  to  sufficient 
accuracy.  Sin'  '  the  method  is  derived  using  the  variational  principle,  It  is  a  second 
order  method  and  converges  rapidly. 


3.2.2  The  initial  guess  for  the  eigenvalue 

The  second  order  eigenvalue  Iteration  procedure  described  above  will  only 
work  if  the  trial  mode  function  has  the  correct  number  of  zero  crossings.  Thus, 
some  other  method  must  be  used  to  get  an  approximate  value  for  the  eigenvalue  in 
the  correct  range.  For  the  first  mode  the  WKB  method  could  be  used  to  get  an 
approximate  value  for  kQ,  or  a  crude  trial  function  could  be  used  with  the 

variational  principle.  An  even  cruder  method  is  simply  to  take  u/c_  and  o/c  as 

o  m 

limits  for  k  and  to  use  the  shooting  method  together  with  a  binary  search.  If 
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there  are  too  many  zero  crossings  of  the  mode,  the  trial  value  of  k  is  too  small;  if 

n 

there  are  too  few  zero  crossings,  the  trial  value  is  too  large.  The  range  of  possible 
values  for  k^  is  halved  at  each  iteration.  Once  the  number  of  mode  crossings  is 

correct  the  faster  second  order  method  can  be  used. 

For  higher  modes  the  upper  limit  on  k^  is  the  wave  number  of  the  mode 

(n-1).  The  lower  limit  is  not  as  easy  to  determine,' although  q/c_  could  always  be 

23 

used.  For  an  isovelocity  profile,  which  is  the  worst  case,  an  approximation  to  k  is 

n 

given  by: 

0  2 

k2  =  k  2  -  _  (n-2)  .  (29) 

n  n-i  .  2 


Thus,  if  we  take  the  upper  limit  from  the  above  equation,  with  a  safety  factor  of 

two  for  good  measure  multiplying  the  second  term  on  the  right,  we  should  be  safe. 

A  conservative  estimate  for  the  lower  limit  of  k  is  thus: 

n 


kmln 

n 


max  [k  2  - 
CB  0-1 


~2(n-2)]  l/t } . 
h2 


(30) 


The  trial  value  of  k^  for  higher  modes  can  be  an  extrapolation  procedure 

based  on  the  values  of  k  for  the  previous  two  or  three  modes. 

n 


3.2.3  The  number  of  trapped  modes 

At  a  given  frequency  the  number  of  modes  trapped  by  a  given  profile  can  be 
determined  by  applying  the  shooting  method  once,  with  a  trial  value  of  k^  equal  to 

o/Cg.  The  number  of  zero  crossings  are  then  counted  and  the  correction  of  Eq. 

(28)  applied  to  k^.  If  the  new  value  of  k^  is  greater  than  u/Cg  then  the  mode  is 

trapped;  i.e.  there  are  N  modes.  If  the  new  value  of  k^j  is  less  than  o/Cg,  only  N-1 

modes  are  trapped. 

This  procedure  can  also  be  used  to  determine  which  modes  have  a  phase 
velocity  less  than  some  specified  phase  velocity.  By  applying  the  procedure  twice, 
the  number  of  modes  between  two  phase  velocities  can  be  calculated. 
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3.2.4  Prevention  of  exponent  overflow 

In  the  shooting  procedure,  the  integration  limits  z  and  z.  must  be  chosen  so 

ft  D 

that  the  function  u  (z)  is  effectively  zero  for  z  <  z  and  z  >  zt .  This  can  be  done 
n  a  b 

approximately  by  evaluating  the  WKB  integral: 

ICz)  =  /  2  |y(z)|  dz  (31) 

z. 


where  z  is  one  of  the  turning  points  z  or  z  ,  y  bas  been  defined  earlier,  and  the 

l  12 

integral  is  evaluated  upward  from  z  or  downward  from  z  .  The  numerical 

1  2 

integration  does  not  have  to  be  performed  with  any  great  accuracy,  so  the 
trapezoidal  rule  is  sufficiently  accurate.  When  the  value  of  the  integral  has 
reached  some  value,  say  D,  then  the  mode  function  will  be  reduced  in  amplitude  by 

a  factor  of  exp(D)  compared  with  its  value  at  the  turning  point.  The  depth  z  or  z. 

a  b 

at  which  the  integral  exceeds  D  is  then  used  for  the  starting  point  for  the 
(accurate)  numerical  integration  of  the  differential  equation. 


3.2.5  Initial  values  for  the  shooting  procedure 

To  obtain  mode  functions  which  have  an  approximate  value  of  unity  at  the 

turning  points,  suitable  values  of  u  and  u'  should  be  chosen  at  z  and  z. . 

a  d 

Near  the  surface  or  z  : 

a 

u(z  )  =  mint - - - ,  1}  if  z  ■  0  and  z  >  0  (32) 

«•  sinh[I  (z )]  *  i 

i  a 


u(z  )  =  expt-I  (z  )]  if  z  >  0  (33) 

ft  X  ft  ft 

In  either  case, 

u'(z  )  »  y(z  )  u(z  ).  (34) 

ft  ft  ft 

In  the  bottom,  since  the  solution  behaves  as  exp[-I  (z)],  the  starting 
conditions  are: 


ufz^)  *  expl-Hz^)] 

(35) 

u'fZj)  »  -ylz^)  u(zfe)  ‘ 

(36) 

J 


■j 

'1 
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3.3  Discussion 

If  the  speed  of  sound  at  the  surface  is  greater  than  the  phase  velocity  of  the 

mode;  i.e.  [u/c(o)]<k  ,  then  the  two-ended  shooting  method  is  superior  to 

n 

one-ended  shooting  methods.  Alsof  if  the  sound  speed  profile  is  such  that 

[o/c(0)]>k  ,  then  the  two-ended  shooting  method  offers  no  disadvantages  over  the 
n 

one-ended  shooting  methods;  (unless,  the  eigenvalue  refinement  procedure  of  Eq. 
(28)  happens  to  be  superior).  If  the  sound  speed  profile  has  two  or  more  local 
minima,  as  would  occur  for  example  in  a  surface  duct,  then  for  modes  which  have 
several  turning  points,  loss  of  precision  will  occur  in  integrating  between  the  two 
minima.  This  loss  of  precision  will  be  about  the  same  no  matter  how  we  integrate 
through  the  sound  speed  maximum,  whether  by  a  one-  or  a  two-ended  shooting 
method.  However,  the  two-ended  method  has  the  advantage  that  any  loss  of 
precision  near  the  surface  is  avoided,  and  moreover  the  method  described  in  Sec. 
3.2  has  promise  for  isolating  the  two  regions  if  necessary.  Thus,  there  seems  to  be 
no  penalty  in  using  the  two- ended  shooting  method  and  for  many  sound  speed 
profiles  greater  precision  and  stability  can  be  achieved  compared  to  the 
conventional  single-ended  shooting  methods. 


4.  RESULTS 

Two  computer  codes  have  been  developed  which  use  the  two-ended  shooting 
technique.  One  PROLOS/MODES  [EUis  and  Leverman,  1982;  Leverman,  1982]  uses 
a  layered  environmental  model  similar  to  Bartberger  [1973];  PROLOS  has  been  used 
extensively  at  DREA  since  1979.  A  more  recent  program  PRLOSS/NORMOD 
[MacBachem,  1983]  uses  a  combination  of  the  Noumerov  method  and  a 
Runge-Kutta  ordinary  differential  equation  solver. 

In  this  section  we  present  some  results  based  on  using  these  models:  the 
AESD  Workshop  Test  Cases,  mode  loss  calculations,  group  velocity  calculations, 
and  comparison  with  shallow  water  propagation  loss  measurements. 


4.1  AESD  Workshop  Test  Cases 

Figs.  3a-3c,  show  the  profiles  used  as  test  cases  for  the  AESD  Workshop  on 
low  frequency  modelling  [Spofford,  1973].  Some  of  them  are  quite  difficult  tests  of 
normal  mode  programs.  The  original  tests  had  units  of  feet  and  feet  per  second; 
here  they  have  been  converted  to  metric  using  the  exact  conversion  factor  of 
0.3048  m  per  foot. 
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Fig.  3.  Sound  speed  profiles  for  the  AESD  Workshop  Test  Cases. 


AESD  TEST  CASE  1 
PACIFIC  PROFILE 


3(a)  -  Test  Case  1:  Pacific  profile. 


AESD  TEST  CASE  2 
NORTHEAST  ATLANTIC  PROFILE 


3(b)  -  Test  Case  2:  Northeastern  Atlantic  profile. 


Fig.  3  (Cont'd) 


AESD  TEST  CASE  3 
SHALLOW  WATER 
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Sound  Speed  (m/s) 


3(c)  -  Test  Case  3:  Shallow  water  profile. 


Figs.  4a-4e  show  the  propagation  loss  calculations  for  five  test  cases:  two 
frequencies  300  and  25  Hz  for  Test  Case  1,  one  frequency  for  Test  Case  2,  and  two 
frequencies  500  and  50  Hz  for  Test  Case  3.  Calculations  were  performed  using  the 
PROLOS/MODES  normal  mode  programi  with  500  homogeneous  layers  of  equal 
thickness  in  the  water  column.  The  sound  speed  used  in  each  layer  was  obtained  by 

—  2 

first  interpolating  between  the  input  profile  points  so  that  c  (z)  was  linear,  then 
using  the  average  value  of  the  interpolated  sound  speed  in  each  layer.  No  earth 
curvature  correction,  absorption  losses,  or  other  losses  are  included  in  the 
calculations.  Calculation  times  were  about  4  to  5  seconds  per  mode  on  a 
DEC-2060  computer  for  eigenvalue  accuracies  of  15  digits.  Typically  only  about  5 
iterations  of  the  second  order  metnod,  Eq.  (28),  are  required. 

As  mentioned  above,  some  of  these  test  cases  are  quite  difficult  tests  of 
normal  mode  models,  and  there  is  considerable  variation  in  the  results  presented  at 
the  AESD  Workshop  [Spofford,  1973].  Indeed  many  of  the  normal  mode  programs 
could  not  handle  all  the  tests.  There  do  not  seem  to  be  benchmark  solutions  for 
most  of  these  cases  but  our  predictions  are  in  reasonable  agreement  with  other 
models.  The  only  difficulty  we  encountered  with  these  profiles  was  in  Test  Case  2, 
the  northeast  Atlantic  profile,  which  has  two  major  sound  channels  separated  by  a 
local  sound  speed  maximum.  In  this  case  the  mode  function  4  fails  to  converge 
after  30  Iterations  although  the  wave  number  seems  accurate.  The  problem  seems 
to  be  due  to  the  fact  that  the  mode  has  a  large  amplitude  in  the  upper  duct,  and  a 
small  amplitude  in  the  lower  duct  where  the  matching  depth  was  chosen. 


4(a)  -  Test  Case  1A:  frequency  300  Hz,  source  depth  91.44  m,  receiver  depth 
27.432  m.  The  solid  line  is  for  the  surface  channel  only  with  a 
half-space  below;  the  dashed  line  is  for  the  entire  profile  but  with 
modal  phase  velocities  confined  to  the  range  1536  to  1539.24  m/s. 


PACIFIC  PROFILE  -  25  Hz 


Fig.  4  (Cont’d) 


NORTHEAST  ATLANTIC  PROFILE  -  50  Hz 
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4(c)  -  Test  Case  2:  frequency  SO  Hz,  source  depth  243.84  m, 
receiver  depth  1097.28  m. 
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4(d)  -  Test  Case  3A:  frequency  500  Hz,  source  depth  6.096  m, 
receiver  depth  12.192  m.  The  dashed  line  shows  the  loss 
when  the  inodes  are  summed  incoherently. 


Fig.  4  (Cont'd) 


SHALLOW  WATER  -  50  Hr 


4(e)  -  Test  Case  3B:  frequency  SO  Hz,  source  depth  6.096  m, 
receiver  depth  12.192  m. 

Of  particular  interest  is  Case  1A,  the  Pacific  surface  duct  at  300  Hz,  where 
normal  mode  methods  would  be  expected  to  have  difficulties.  Calculations  were 
performed  using  the  surface  duct  only  with  a  half-space  (p  »  1,  c  »  1539.24  m/s) 
below  152.4  m,  as  well  as  for  the  entire  profile  with  phase  velocities  confined  to 
the  range  1536  m/s  to  1539.24  m/s.  Fig.  4(a)  shows  that  the  effect  of  using  the 
whole  profile  is  to  reduce  the  propagation  loss  somewhat,  and  to  cause  the  nulls  to 
be  somewhat  filled  in.  No  computational  difficulties  were  encountered  in  spite  of 
the  two-duct  nature  of  the  profile  and  the  higher  frequency  than  Case  2. 

There  is  a  benchmark  result  [Stickler,  1975]  for  Case  3B  where  there  is  only 
one  trapped  mode,  near  cutoff.  In  this  case,  the  continuous  spectrum,  neglected  in 
our  model,  is  important.  Our  results  agree  with  Stickler's  results  when  he  excludes 
the  contribution  of  the  continuous  spectrum. 


4.2  Mode  Attenuation  Coefficients 

The  AESD  test  cases  discussed  in  Sec.  4.1  did  not  include  any  loss 
mechanisms,  other  than  geometrical  spreading.  Fig.  5  shows  the  attenuation 
coefficients  of  mode  1  as  a  function  of  frequency  for  the  AESD  Test  Case  3.  The 
following  parameters  were  used  for  the  calculations: 
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Bottom  absorption  coefficient  a  =0.2  dB/m-kHz 
Surface  roughness  a  =  0.2  m 

Bottom  roughness  a  =  0.2  m 

i 

Shear  speed  in  bottom  c  =  200  m/s 

s 

Shear  wave  absorption  coefficient  a$  =  1  dB/m-kHz. 

The  PRLOSS/ NORMOD  model  [MacEachern,  1983}  was  used  for  these  calculations* 
since  the  individual  modal  attenuation  coefficients  were  readily  available  and  it 
was  possible  to  do  multiple  frequencies  conveniently. 


£ 


Fig.  5.  Attenuation  coefficients  of  mode  1  for  AESD  Test  Case  3  with 
various  assumed  loss  mechanisms:  BOTTOM  CS  ■  0  [Eq.  (13) 
with  cs  =  0];  BOTTOM  CS  =  200  [Eq.  (19)];  WATER  [Bq.  (15)]; 

B-SCATT  [Eq.  (17)];  and  S-SCATT  [Bq.  (16)]. 


Note  that  all  the  losses,  except  the  bottom  loss  calculated  by  the 
perturbation  formula,  Eq.  (13),  approach  0  as  the  frequency  approaches  cutoff  from 
above  50  Hz.  The  surface  scattering  loss  increases  to  a  maximum  near  80  Hz,  then 
decreases  to  0  near  180  Hz  when  the  phase  velocity  of  mode  1  becomes  less  than 
the  sound  speed  at  the  surface.  At  high  frequencies  the  volume  absorption  in  the 

water  and  the  bottom  scattering  loss  behave  as  f 2.  The  bottom  loss  decreases  with 
increasing  frequency  since  the  grazing  angle  becomes  smaller  and  the  mode 
function  does  not  penetrate  as  deeply  into  the  bottom. 
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Notice  that  the  losses  with  shear  waves  included  are  less  than  the  losses 

with  c  =0.  Further  calculations  show  that  the  effect  of  shear  waves  is  negligible 
s 

here  since  the  bottom  is  so  lossy  anyway.  The  problem  is  that  the  two  formulae  for 
calculating  bottom  losses,  Eqs.  (13)  and  (19),  are  inconsistent,  a  problem  that  needs 
further  investigation. 


4.3  Group  Velocities  and  Mode  Functions 

The  group  velocity  rather  than  the  phase  velocity  determines  the  speed  at 
which  energy  is  transported  for  a  broadband  signal.  For  a  multilayered 
environment,  the  group  velocity  of  each  mode  can  show  a  rather  complicated 
frequency  dependence  with  several  local  minima  and  maxima.  Chapman  and  Ellis 
[1983a]  have  shown  that  the  mode  functions  can  be  used  to  give  an  intuitive 
understanding  of  the  group  velocity. 

Figs.  6a  and  6b  [Chapman  and  Ellis,  1983a]  show  group  velocity  curves  and 
mode  functions  for  mode  5  in  a  three  layer  model:  water,  p  =  1,  h  *  100  m,  c  ■ 

ill 

1500  m/s;  silt,  p  =  1,  h  »  100  m,  c  *  1650  m/s;  rock,  p  =  2,  h  =  «,  c  =  4500 

2  2  2  9  2  2 

m/s,  where  the  h^  refer  to  layer  thicknesses.  The  maxima  in  the  group  velocity 

curves  at  50  and  73  Hz  correspond  to  a  "resonance"  in  the  silt  layer.  The  maximum 
near  cutoff  at  18.92  Hz  occurs  because  most  of  the  normalization  of  the  mode 
function  occurs  in  the  high-speed  rock  layer.  The  minima  near  20  and  90  Hz  occur 
when  the  equivalent  ray  angle  in  the  silt  or  water  is  nearing  the  critical  angle  with 
the  higher  speed  layer  below.  The  higher  speed  layer  does  not  contribute  enough  to 
increase  the  group  velocity  however. 

The  group  velocity  and  mode  functions  were  calculated  using  the 
PR OLO S/MODES  program  and  some  small  utility  programs. 


Fig.  6a.  Group  velocity  versus  frequency  for  mode  5  for  a  three-layer 
model.  Symbols  refer  to  mode  plots  in  Fig.  6b. 
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Fig.  6b.  Amplitude  of  mode  5  as  a  function  of  depth  for  the  three-layer  model. 
The  frequencies  correspond  to  the  points  marked  (x)  in  Fig.  6a. 


4.4  Comparison  With  Shallow  Water  Data 

It  is  usually  much  more  difficult  to  obtain  good  agreement  with 
experimental  results  than  with  other  models.  In  shallow  water  the  main  difficulty 
is  in  obtaining  accurate  environmental  information,  particularly  about  the  sea  bed, 
for  input  to  the  model. 

Fig.  7  shows  a  comparison  of  propagation  loss  predictions  with 
measurements  at  a  shallow  water  site.  The  measurements  were  made  in  summer 
conditions  in  a  range-dependent  environment  with  water  depth  varying  between 
115  and  155  m,  and  a  sandy  sediment  layer  overlying  a  semi-consolidated  bottom. 
The  sediment  properties  were  obtained  from  an  independent  geo-acoustic  model 
based  on  the  analysis  of  a  vertical  incidence  seismic  profiler.  The  calculations 
were  done  by  the  adiabatic  range-dependent  normal  mode  model  PROLOS.  There 
is  good  agreement  between  the  model  and  measurements  at  practically  all  ranges 
and  frequencies.  This  shows  about  the  best  agreement  that  one  is  likely  to  obtain 
in  shallow  water,  and  is  especially  gratifying  since  it  was  not  necessary  to  tune  the 
geo-acoustic  parameters  to  improve  the  fit.  More  details  about  the  model  and 
experiment  are  given  by  Ellis  and  Chapman  [1984]. 


RANGE  (km) 

Fig.  7.  Measured  and  modelled  propagation  loss  versus  range  for  a  range 
dependent  shallow  water  site.  Measured:  □  64  Hz,  A  256  Hz, 

0  1024  Hz.  Modelled:  solid  lines. 


4.5  Discussion 

In  addition  to  the  results  shown  earlier  in  this  section,  the  two-ended 
shooting  technique  has  been  used  extensively  in  shallow  water  calculations  at 
DREA  since  1979;  in  addition  to  the  papers  quoted  earlier  in  this  section  see:  [Ellis 
and  Chapman,  1980],  [Chapman,  1980],  [Chapman,  1983],  [Chapman  and  Ellis,  1983] 
as  well  as  various  other  DREA  reports.  The  method  has  worked  successfully  on  all 
shallow  water  profiles  and  bottom  models  that  we  have  tried  and  at  frequencies  up 
to  3000  Hz.  This  is  not  to  say  that  we  have  always  been  able  to  get  good 
agreement  with  our  experimental  measurements;  but  rather  that  the  environmental 
inputs  are  inaccurate  or  the  model  does  not  account  for  certain  phenomena. 

One  specific  thing  that  is  lacking  is  the  ability  to  handle  shear  waves  when 
the  bottom  is  hard  (e.g.  granite  with  c  ~  3000  m/s)  or  consolidated  (e.g.  limestone 

or  chalk  with  c&  ~  1000  m/s).  In  the  former  case  our  algorithm  could  easily  be 

generalized  to  handle  a  bottom  half -space  since  it  would  only  involve  a  change  in 
the  boundary  condition  at  z^.  In  the  latter  case  the  shear  speed  is  too  high  for 

perturbation  techniques  to  be  applied,  and  the  mode  functions  and  wave  numbers 
would  be  complex  since  c$  <  c(z)  for  z  in  the  water.  See  [Ellis  and  Chapman,  1984] 

for  more  discussion. 
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In  deep  water  the  model  neglects  the  continuous  spectrum  (or  equivalently 
the  steep  rays)  which  may  be  important  at  short  ranges.  As  well,  the  computation 
time  increases  since  the  number  of  modes  is  proportional  to  water  depth.  Other 
techniques  such  as  the  fast  field  method  or  parabolic  equation  may  be  more 
attractive.  However,  if  long  range  propagation  results  are  required,  normal  mode 
theory  is  often  adequate  or  preferable.  In  such  cases  the  algorithm  would  be  quite 
useful,  especially  if  a  restricted  range  of  phase  velocities  is  all  that  needs  to  be 
computed.  Since  a  two-ended  shooting  technique  is  used,  the  high  sound  speed  near 
the  surface  is  not  a  problem,  as  would  be  the  case  in  single-ended  shooting 
methods;  the  only  problem  that  seems  to  arise  (in  range-independent  cases)  is  in 
the  two  duct  problems  such  as  AESD  Case  2  where  there  are  two  strong,  well 
isolated  sound  channels.  In  such  cases  a  multiple  shooting  technique  or  a  finite 
difference  approach  might  have  to  be  used  rather  than  the  method  described  here. 


5.  SUMMARY  AND  CONCLUSIONS 

We  have  described  a  two- ended  shooting  method  that  has  been  used  for 
normal  mode  calculations  at  DREA  since  1978.  It  is  a  technique  borrowed  from 
quantum  mechanical  calculations  and  generalized  to  handle  density  changes. 

For  completeness  we  have  included  loss  mechanisms  needed  to  compare 
model  predictions  with  experiment,  and  have  shown  a  variety  of  results.  These 
include  propagation  loss  calculations  for  standard  test  cases,  mode  functions  and 
group  velocities,  model-data  comparisons,  and  attenuation  losses  due  to  surface 
and  bottom  volume  absorption,  surface  and  bottom  roughness,  and  shear  wave 
losses  in  the  bottom.  Ways  in  which  the  method  could  be  extended  or  improved  are 
also  discussed. 

The  method  has  been  successfully  applied  to  all  shallow  water  profiles  that 
we  have  attempted,  and  at  frequencies  up  to  3000  Hz.  We  have  not  always 
obtained  good  agreement  with  experimental  measurements,  but  this  seems  to  be 
due  to  inaccuracies  in  the  environmental  information  or  the  fact  that  our  model 
does  not  account  for  certain  things  (e.g.  high  speed  shear  waves  in  the  bottom).  i 

I 

i 

In  deep  water  environments  the  model  usually  performs  well,  although  it  of 
course  neglects  the  continuous  spectrum  (or  equivalently,  the  steep  rays)  which 
may  be  important  at  short  ranges.  Well  isolated  sound  channels  (such  as  the  AESD 
Workshop  Test  Case  2)  can  cause  some  of  the  mode  functions  to  be  inaccurate  at  j 

frequencies  above  100  Hz.  On  the  other  hand,  the  algorithm  using  the  layered  i 

environmental  model  has  no  trouble  with  AESD  Test  Case  1A,  a  surface  duct  at  300  i 

Hz,  even  when  the  complete  profile  is  used. 
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Appendix  A:  Notation 
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cm 

CB 
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PL 

IPL 


space  co-ordinate  (x,y,z)  or  (r,z,0) 

position  of  source  (0,z  ,0) 

o 

range  co-ordinate 
depth  co-ordinate 
density 

density  of  water 
sound  speed 

frequency,  angular  frequency  o>*2uf 

time 

pressure 

mode  function 

water  depth,  layer  thickness 

Hankel  function  of  zeroth  order  and  first  kind 
wave  number 

wave  number  of  the  n-th  discrete  normal  mode 

imaginary  part  of  the  wave  number  of  the  n-th  mode 

complex  wave  number  of  the  nth  mode;  K  *  k  +  16 

n  n  n 

minimum  sound  speed 

sound  speed  in  bottommost  layer 

source  strength 
propagation  loss  in  dB 
incoherent  propagation  loss  in  dB 
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Kronecker  delta 


6 

mn 

a  attenuation  coefficient  (dB/length  at  1  kHz) 

c  attenuation  coefficient  (nepers  Aength) 

a  ,o  r.m.s.  heights  of  rough  surface  and  bottom 

o  1 

R  plane  wave  reflection  coefficient 

Y  (z)  local  vertical  wave  number 


0  .  0  , 

1  2  S 

(n) 

vg 

I(z) 

d> 

z  ,z 
1  2 


vertical  wave  numbers  used  to  determine  reflection  coefficient 

group  velocity  of  mode  n 

phase  or  exponent  integral 

phase  shift  at  a  mode  turning  point 

upper  and  lower  turning  points  of  mode  function 


z  ,z.  upper  and  lower  limits  for  the  numerical  integration  of  the  d.e. 

a  b 


z  matching  depth  for  the  two-ended  shooting  method 

m 


Dimensionless  Variables 


x  »  z/h 


V(x)  =  [hu/c(z)] 


_  .2,  2 
Q  3  b  k 
n  n 


.v, 

w  *  h  u 
n  n 

dw  .  du 

w  =_a  =h,/2  — - 
n  dx  dz 


s(x)  «  p(z)/p 


normalized  depth 
dimensionless  sound  speed 
dimensionless  eigenvalue 

dimensionless  eigenfunction 

dimensionless  derivative 
dimensionless  density 
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Appendix  B:  Absorption  as  a  Perturbation 


If  absorption  is  introduced  into  the  medium,  the  plane  wave  becomes 
attenuated;  that  is, 


ikr  ikr  -er 
e  ->  e  e 


(B-l) 


where  e  is  the  attenuation  coefficient  in  nepers  per  unit  distance.  This  is 
equivalent  to  letting  the  local  wave  number  of  the  sound  speed  become  complex; 
that  is, 

k(z)  =  o/c(z)  +  i  e(z)  .  (B-2) 

The  corresponding  normal  mode  wave  number  also  becomes  complex: 

K  =  k  +16.  (B-3) 

n  n  n 

The  imaginary  part  of  6  can  be  obtained  by  a  perturbation  calculation,  in 

n 

which  the  squares  of  small  quantities  are  neglected.  We  start  with  the  normal 
mode  equation: 


du  (z) 

p(z)  A_[  - 5 — ]  +  [k2(z)  -  K  2)  u  (z)  -  0  . 

dz  p(z)  dz  n  n 


CB-4) 


First,  we  multiply  Eq.  (B-4)  by  [p(z>]  un*(z),  where  uq*  is  the  complex  conjugate  of 

ufl.  Then  we  take  the  complex  conjugate  of  Eq.  (B-4)  and  multiply  it  by  [p(z)]  ^^(z). 
Next,  subtract  the  two  equations  and  integrate  from  zero  to  infinity,  giving: 

.  .  du  j  i  du*(z) 

J°°tu*(z)_i  (-i - ?)  -  u  (z)-i  ( J - —  )]  dz 

°  n  dz  p(z)  dz  n  dz  p(z)  dz 


,  .  u  (z)u*(z) 

+  S  “«ka(z)  -  k*\z)]  -  (K  2  -  K  *  )}  — - dz  .  0. 

o  n  n  ni’t) 


The  integrand  of  the  first  equation  can  be  written: 


(B-5) 


,  u*(z)  du  u  (z)  du*(z) 

_d  ^  _n _  _ n  _  n 

dz  p(z)  dz  p(z)  dz 


(B— 6) 


and  since  u  (0)*0,  and  lim  u  (z)  -*  0,  the  integral  is  zero.  Also,  if  the  imaginary 
n  n 

parts  of  k(z)  and  K  are  small, 
n 


v>'  -  V 


CB-7) 


K  2  -  K*  *=  4i  k 
n  n  n 


6 

n 


k2(z)  -  k*  (z)  =  4i  oc( z)/c(z)  . 


Therefore  from  Eqs.  (B-5MB-7),  we  obtain: 


e(z)u  (z)  u  (z) 

f°°  _ 5 _ - —  dz 

£  _  0  p(z)  c(z) _ 

n  k  u  (z)  u  *(z) 

n  p«>  n  n 


p(z) 


dz 


(B-8) 


(B-9) 


If  the  mode  functions  are  approximated  by  their  unperturbed  values,  the 
integral  in  the  denominator  is  Just  the  normalization  (unity).  Moreover,  if  the 
unperturbed  mode  functions  are  used,  the  integrand  of  the  first  integral  in  Eq. 
(B-5)  is  zero  everywhere,  and  not  only  at  the  end  points.  The  attenuation  due  to 
absorption  in  any  layer  is  given  by 


6l.-2-X  1 

"  kn  *1-1 


e(z)  u2(z) 

_ ^_dz. 

p(z)  c(z) 


(B-10) 


For  bottom  absorption  z^  ^  =  h,  and  z^  *  «». 

The  above  formula  can  also  be  used  to  calculate  the  volume  absorption  due 
to  the  water,  by  using  0  and  h  as  the  limits  of  the  integration. 

For  an  infinite  half-space  in  the  bottom  (i.e.  the  bottommost  layer), 

u  (z)  =  u  (h)  exp[-y  (z-h)],  z>h  (B-ll) 

n  n  n 


ST  u  a(z)  dz  * 

h  n 


u  2(h) 
n 


c  B 

\^BCB 


(B-12) 


(B-13) 
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Appendix  C:  Iteration  Procedure  for  the  Two-Ended  Shooting  Method 


The  variational  principle  tells  us  that  if  we  have  a  differential  equation  of 
the  form: 


H(x)  u  (x)  =  X.  u  (x) 
n  n  n 


(C-I) 


where  H  is  some  differential  operator,  and  vq  is  a  trial  eigenfunction  satisfying  the 

appropriate  boundary  conditions,  then  an  approximation  Q  to  \  is  given  by: 

n  n 


r  _L  V  (x)  H(x)  V  (x)  dx 


/ - v  (x)  dx 

p(x)  n 

where  p(x)  is  the  appropriate  weighting  function  and  the  integral  is  performed  over 
the  interval  of  the  boundary  value  problem. 

In  the  two-ended  shooting  method,  we  have  a  trial  function  v  (z)  which 

n 

satisfies  the  normal  mode  equation: 


{  —  [  —  — ]  +  —  }  v  (z)  *  Q 

dz  p(z)  dz  p(z)  c**>  n  “  p(z) 


v  (z) 
n 


(C-3) 


on  the  two  intervals  o  S  z  <  z  ,  and  z  >  z  ,  that  is,  everywhere  except  at  z  .  Q 

mm  m  n 

is  the  trial  approximation  to  the  solution  k2.  The  function  v  (z)  is  continuous,  but 

n  n 

[p<z)]  *v  ‘(z)  is  discontinuous  at  z  =  z  .  By  taking  H(x)  to  be  the  quantity  in 
Eq.  (C-3),  we  can  use  Eq.  (C-2)  to  obtain  a  better  approximation  Q 

n 


Q  '  »  J  m  v  (z)  H(z)  v  (z)  dz  +  J 00  v  (z)  H(z)  v  (z)  dz 
non  n  +  n  n 

zm 


+  J  v  (z)  H(z)  v  (z)  dz  . 


(C-4) 


For  notational  simplicity  in  Eq.  (C-4)  the  function  v^(z)  has  been  assumed  to  be 


normalized  to  unity,  so  the  denominator  term  is  absent.  Since  H(z)v  (z)  =  Q  v  (z) 

n  n  n 


on  the  intervals  0  <  z  <  z  and  z  >  z  ,  the  first  two  integrals  can  be  evaluated  to 

m  m 


give: 


I  +  I 
1  2 


=  QnJc 

n  o 


v2(z) 

n 

p(z) 


dz 


(C-5) 


The  final  integral  of  Eq.  (C-4),  although  integrated  over  an  infinitesmally  small 
interval,  gives  a  finite  contribution  due  to  the  presence  of  the  derivatives  in  H(z). 
It  can  be  evaluated  using  integration  by  parts  to  give: 


m 


z 

m 


v  (z)  -1  [J— 
n  dz  p(  z) 


dv 


v  (z) 
_2.]dz-_5_ 
dz  p(z) 


dv  z 
n  ,  m 


dz 


I  -  J 


m 


m 


,  dv 

[—]a  dz  . 
z  p(z)  dz 
m 


(C-6) 


Since  the  integral  term  on  the  r.h.s.  of  Eq.  (C-6)  is  negligible,  the  improved 
estimate  Q  is  given  by: 


v  (z  )  { 
n  m 


dv  (z+) 
m  m 


Q  =  Q  + 

n  n 


p(z+) 

m 


dz 


_1 

p(z 


dv  (z  ) 
m  m 

\  dz 


m 


v  a(z) 
n 

p(z) 


dz 


(C— 7) 


In  Eq.  (C-7),  the  denominator  has  been  re-introduced  so  the  expression  can  be  used 
with  trial  functions  which  are  not  normalized. 


Appendix  D:  Dimensionless  Co-ordinates 


It  is  often  useful  for  numerical  purposes  to  define  the  variables  in  terms  of 
dimensionless  quantities.  The  calculations  can  then  be  done  independently  of  the 
units  used. 


A  dimensionless  depth  co-ordinate  can  be  defined  in  terms  of  the  water 
depth  h  as: 

x  =  z/h  .  CD-I) 

If  the  following  definitions  are  made: 

V(x)  =  [ho/c(hx)] 2  (D-2) 


_  ,  2,  2 
Q  =  h  k 
n  n 


s(x)  =  p(hx)/p 


w 


then  the  normal  mode  Eq.  (6)  can  be  written: 


dw 

s(x)  —  [  -i - -)  +  rv(x)  -  Q  ]  w  (x)  »  0 

dx  s(x)  dx  n  n 


(D-3) 

(D-4) 


(D-S) 


where  the  w  are  the  new  normal  mode  functions  normalized  so  that: 
n 

J°°s  ^x)  w  2(x)  dx  =  1  . 
o  n 

In  terms  of  the  dimensionless  mode  functions: 

u  (hx)  *  p~l  h"  ^  w  (x) 
n  w  n 

du  (hx)  .  dw  (x) 

n  —  i  v— a/2  n 

-  =  p  n  -  . 

dz  W  dx 


CD-6) 


CD- 7) 

(D-8) 
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The  constraints  on  the  eigenvalues  become: 


ho  ^  ^  ^  >  q  ^  ho> 

V  1  *  N  CT' 


The  number  of  modes,  of  course,  does  not  change. 


The  boundary  conditions  do  not  change,  so  vMO)  =  0,  and  p 


continuous,  and  the  behaviour  of  w  for  large  x  is: 

n 


w  (x)  ~  exp(-g  x) 
n  n 


where, 


g  .|Q 
n  n  a 

c  B 


The  iteration  formula  for  the  eigenvalues  can  be  written  as: 


w  (x  )  w  (x  ) 


w  (x  )[  *L3-  -  -UL] 

n  m  .  + 


Q  new  _  old  + 

r»  x  ft 


s(x .  )  s(x  ) 
m  m 


(D-9) 


w  ,  are 
n 

(D-10) 


(D-ll) 


(D-12) 
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